A core motivation of science is to evaluate which scientific model best explains observed data. Bayesian model comparison provides a principled statistical approach to comparing scientific models and has found widespread application within cosmology and astrophysics. Calculating the Bayesian evidence is computationally challenging, especially as we continue to explore increasingly more complex models. The Savage-Dickey density ratio (SDDR) provides a method to calculate the Bayes factor (evidence ratio) between two nested models using only posterior samples from the super model. The SDDR requires the calculation of a normalised marginal distribution over the extra parameters of the super model, which has typically been performed using classical density estimators, such as histograms. Classical density estimators, however, can struggle to scale to high-dimensional settings. We introduce a neural SDDR approach using normalizing flows that can scale to settings where the super model contains a large number of extra parameters. We demonstrate the effectiveness of this neural SDDR methodology applied to both toy and realistic cosmological examples. For a field-level inference setting, we show that Bayes factors computed for a Bayesian hierarchical model (BHM) and simulation-based inference (SBI) approach are consistent, providing further validation that SBI extracts as much cosmological information from the field as the BHM approach. The SDDR estimator with normalizing flows is implemented in the open-source harmonic Python package.

We advocate for a new paradigm of cosmological likelihood-based inference, leveraging recent developments in machine learning and its underlying technology, to accelerate Bayesian inference in high-dimensional settings. Specifically, we combine (i) emulation, where a machine learning model is trained to mimic cosmological observables, e.g. CosmoPower-JAX; (ii) differentiable and probabilistic programming, e.g. JAX and NumPyro, respectively; (iii) scalable Markov chain Monte Carlo (MCMC) sampling techniques that exploit gradients, e.g. Hamiltonian Monte Carlo; and (iv) decoupled and scalable Bayesian model selection techniques that compute the Bayesian evidence purely from posterior samples, e.g. the learned harmonic mean implemented in harmonic. This paradigm allows us to carry out a complete Bayesian analysis, including both parameter estimation and model selection, in a fraction of the time of traditional approaches. First, we demonstrate the application of this paradigm on a simulated cosmic shear analysis for a Stage IV survey in 37- and 39-dimensional parameter spaces, comparing $\Lambda$CDM and a dynamical dark energy model ($w_0w_a$CDM). We recover posterior contours and evidence estimates that are in excellent agreement with those computed by the traditional nested sampling approach while reducing the computational cost from 8 months on 48 CPU cores to 2 days on 12 GPUs. Second, we consider a joint analysis between three simulated next-generation surveys, each performing a 3x2pt analysis, resulting in 157- and 159-dimensional parameter spaces. Standard nested sampling techniques are simply not feasible in this high-dimensional setting, requiring a projected 12 years of compute time on 48 CPU cores; on the other hand, the proposed approach only requires 8 days of compute time on 24 GPUs. All packages used in our analyses are publicly available.

Spurio Mancini et al. (2022) (SM22 hereafter), the exploration of high-dimensional parameter spaces in particular, developed CosmoPower, a suite of - O(100) parameters and higher - necessary to accurately neural network emulators of cosmological power spectra model the physical signals and their several systematic that replaces the computation of these quantities traditionally contaminants. Sampling the posterior distribution performed with Einstein-Boltzmann solvers such in these high-dimensional spaces represents a significant as the Code for Anisotropies in the Microwave Background computational challenge for Markov Chain Monte Carlo (CAMB, Lewis & Challinor 2011) or the Cosmic (MCMC) algorithms (Roberts et al. 1997; Katafygiotis Linear Anisotropy Solving System (CLASS, Blas et al. & Zuev 2008; Liu 2009), which are traditionally used 2011). In SM22 the authors show how Bayesian inference in cosmological analyses (Lewis & Bridle 2002; Audren of cosmological parameters can be accelerated by several et al. 2013; Brinckmann & Lesgourgues 2019; Torrado & orders of magnitude using CosmoPower; the speed-up becomes Lewis 2021). Gradient-based inference methods, such as particularly relevant when the emulators are employed Hamiltonian Monte Carlo (HMC, Duane et al. 1987; Neal within an inference pipeline that can be run on 1996) and variational inference (VI, Hoffman et al. 2013; graphics processing units (GPUs). Blei et al. 2017), manage to concentrate the sampling An additional advantage in using machine learning in regions of high posterior mass, even in large parameter emulators is that they efficiently provide accurate derivatives spaces, provided one has efficient access to accurate with respect to their input parameters. This is derivatives of the likelihood function with respect to made possible by the automatic differentiation features the model parameters (Brooks et al. 2011; Neal 2011; implemented in the libraries routinely used to build these Zhang et al. 2017; Betancourt 2017).

Bayesian inference applied to microseismic activity monitoring allows the accurate location of microseismic events from recorded seismograms and the estimation of the associated uncertainties. However, the forward modelling of these microseismic events, which is necessary to perform Bayesian source inversion, can be prohibitively expensive in terms of computational resources. A viable solution is to train a surrogate model based on machine learning techniques, to emulate the forward model and thus accelerate Bayesian inference. In this paper, we substantially enhance previous work, which considered only sources with isotropic moment tensors. We train a machine learning algorithm on the power spectrum of the recorded pressure wave and show that the trained emulator allows complete and fast event locations for $\textit{any}$ source mechanism. Moreover, we show that our approach is computationally inexpensive, as it can be run in less than 1 hour on a commercial laptop, while yielding accurate results using less than $10^4$ training seismograms. We additionally demonstrate how the trained emulators can be used to identify the source mechanism through the estimation of the Bayesian evidence. Finally, we demonstrate that our approach is robust to real noise as measured in field data. This work lays the foundations for efficient, accurate future joint determinations of event location and moment tensor, and associated uncertainties, which are ultimately key for accurately characterising human-induced and natural earthquakes, and for enhanced quantitative seismic hazard assessments.